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1. Introductory remarks and computational setup 

The European Twisted Mass Collaboration (ETMC) has recently performed simulations em- 
ploying mass degenerate light (up/down) doublet quarks and a heavy mass non-degenerate pair 
for strange and charm quarks entering in the era of precise and realistic Nf ■ = 2 + 1 + 1 dynam- 
ical quark lattice computations [JIJ]. In the ETMC setup gluon interactions are described by the 
Iwasaki action [0]. Fermions are regularised in the maximally twisted mass (Mtm) Wilson lattice 
formulation [Q, [|, |5|]. This choice of the fermionic action has the benefit of achieving automatic 
0(a) -improvement, generally leading to small 0(a 2 ) lattice artefacts. Up to now ETMC has pro- 
duced Nf = 2 + 1 + 1 dynamical quark gauge configurations at three values of the lattice spacing 
(namely a ~ 0.06,0.08 and 0.09 fm). Gauge ensembles have been produced at several values of the 
quark masses with the lowest pseudoscalar mass being of about 270 MeV. Obviously the inclusion 
of the dynamical strange and charm quarks offers the advantage of taking into account a rather 
important source of systematic effects. 

Computation of renormalisation constants (RCs) is a crucial step in order to extract physical 
quantities from lattice data. It is worth noticing that the Mtm setup offers a rather convenient quark 
mass renormalisation pattern. Indeed, for the renormalisation of the masses of degenerate and non- 
degenerate quark pairs one only needs to know the non-singlet Zp and Zj renormalisation constants 
(and not Z s o) M\- RCs of operators with non-zero anomalous dimension need to be computed in 
the chiral limit and for this reason dedicated lattice simulations employing Nf = 4 light and (for 
simplicity) degenerate dynamical quarks are required. Indeed we have produced Nf = 4 gauge 
configuration ensembles corresponding to several sea quark mass values and we determined RCs 
extrapolating their lattice estimators to the chiral limit. We employed the RI-MOM scheme [||] and 
the techniques already used for Nf = 2 RCs [f7|]. However, for the case of Nf = 4 simulations with 
the action and at the lattice spacing values we are currently using, the implementation of maximal 
twist (i.e. tuning the PCAC quark mass to zero), which would guarantee O(a) improvement of RCs, 
is not a trivial task. In fact in the region of small PCAC quark mass values simulation instabilities 
occur that lead for it to very large autocorrelation times. Hence we opted for an alternative way, 
already proposed in Ref. [g], to achieve O(a) improvement though working out of maximal twist. 
The method, being based on averaging results obtained at opposite values of the PCAC quark mass, 
entails the need of doubling the reasonably low CPU time cost for producing gauge simulations 
at non-zero standard and twisted quark mass. The present contribution is a report of a work in 
progress. A first presentation and numerical test of our method appeared already in Ref. [@J. 



1.1 Computational setup 

We consider the following fermionic lattice action (written in the so called physical basis) 



S f = a 4 £ £ q f \y ■ V - i Kr ^(~rV + m cr ) +M, 



x f=l 



0,f 



(l.i) 



where qf denotes a singlet quark flavour and rf takes values either -1 or +1. The chiral quark field 
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rotation Xf — > <lf = ex p[2(f — @o,f)Y5 r f]Xf brings the action into the so-called twisted basis 1 



S" n 



i " fl4 IL^ y-V--V*V+m ,/ + j75W Zf(x) 

X f=l 



(1.2) 



The bare mass parameters and the angle f? are given by 

sin 0o,/ 



M 0,f = V ( m o, 



m. 



cr; 2 +M/ • 



WIt- 



COS do J 



(1.3) 



M ,/ ,J M 0i/ 

In practice we make use of /m PC ac to estimate (otq/— m cr ). In this way we take as the renormalised 
quantities the polar quark mass Mf 



Zp'Mf 



Z P y Z\m^ CAC + nj and the angle df, comple- 
mentary to the twisted angle ft)/ (df = n/2 — (Of), given by tan df = Zawipcac/A 4 /- As we use four 
mass degenerate quark flavours and we adopt a partially quenched setup, the knowledge of the four 
parameters M sea , sea ,M va i, va i is sufficient to describe our RC computation. 

In the following we focus on the evaluation of the RCs of the (non-singlet) quark bilinear 
operators 2 Or = Xf^Xf where F = S,P,V,A,T, in the RF variant of the RI-MOM scheme. One 
first determines the quark field RC, Z q , through 



1 y / 
l2N(p)L 



Tr(Y P Sf(p) 
Pp 



1 . 



any/ 



(1.4) 



where = ±sinap M , p 2 = £^p 2 and S f (p) = a 4 £ A .e ipx (Xf(x)Xf(ty) is the Landau gauge 
quark propagator in momentum space. The sum £ p ' runs over the Lorentz indices for which p p is 
different from zero and N{p) = '1. Then one computes the RC, Zp, of the operator Or via 



7 -i 7 (//') 



Tr 



^(p)G { / f \p,p)S f 1 (p))Pr 



where 



4 ff \p,p)=a^e 



-ip(x-y) 



(Xf(x){x f rxf)(o)Zf(y)) 



S,P,V,A,T. 



(1.5) 



(1.6) 



.v.y 



with Tf = —rf for the Wilson parameters of the (valence) quark flavours / and f . We note that 
RCs are blind to the choice of sign(rf) but lattice artefacts in their estimators in general are not. 

In our computation we will exploit the fact that the 0(a 2k+l ) artefacts occurring in the vac- 
uum expectation values of (multi)local operators O vanish if we take the d -average defined by 



,3 



The 0(a) improvement obtained in this way is a consequence of the 

symmetry 1 3° x [do — > — do) x ^ x (Mo — > —Mo) of the lattice action and occurs for operator 
expectation values and form factors that are invariant under J 2 x (0 O — > — do), see refs. In 
particular this holds for our RCs estimators at any value of quark mass, Mf, and momentum, p. 



'For consistency with the standard Wilson fermion notation, the operator RCs are named according to the form the 
operators take in this basis, where the Wilson term is untwisted. 

2 For the computation of the RCs of the four-fermion operators using the same setup see Ref. [juj. 
3 We denote by £P the parity transformation of the fields and by &j the transformation defined from @^qr(x) = 



e^l 2 q f {-x), %q f (x) = e^l z q f {-x) and 9 d Up = U^-x-afl). 
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2. Analysis and Results 

We have produced Nf = 4 dynamical quark gauge configurations at three values of the inverse 
gauge coupling, j3 = 1.90, 1.95 and 2.10, and for each j3 at a number of values of M sea and (nearly) 
opposite values of sea (the p/m in the ensemble labels refers to sign(6 sea )). An overview of these 
ensembles and the valence mass parameters chosen for the Landau gauge correlation functions is 
given in Table [[} Some quark propagator computations are still in progress, as indicated. Hence, 
the corresponding ensembles can not be used in this analysis. 



CllaCIIlDlC 




" m PCAC 


,,/usea 
CIM Q 


osea 

(7 




"™PCAC 








n 

p = 


1.90 






4m 


0.0080 


-0.0390(01) 


0.0285(01) 


-1.286(01) 


in progress 




4p 


0.0080 


0.0398(01) 


0.0290(01) 


+1.291(01) 


in progress 




3m 


0.0080 


-0.0358(02) 


0.0263(01) 


-1.262(02) 


in progress 




3p 


0.0080 


. n a t r /" / n o \ 

+0.0356(02) 


0.0262(01) 


+1.260(02) 


in progress 




2m 


0.0080 


-0.0318(01) 


0.0237(01) 


-1.226(02) 


in progress 




2 P 


n nnon 

0.0080 


. n no 1 A/A^ \ 

+0.0310(02) 


n noo i /ni \ 

0.0231(01) 


.1 o 1 o /no \ 

+1.218(02) 


in progress 




lm 


0.0080 


-0.0273(02) 


0.0207(01) 


-l. 174(03) 


in progress 




IP 


n nnon 

0.0080 


. n A<in r/A(\ 

+0.0275(04) 


n a^aa /n 1 \ 

0.0209(01) 


. 1 i t~i / nr\ 

+1.177(05) 


in progress 










n 

P = 


1.95 






lm 


0.0085 


-0.0413(02) 


0.0329(01) 


-1.309(01) 


[0.0085,..., 0.0298] 


-0.0216(02) 


IP 


0.0085 


+0.0425(02) 


0.0338(01) 


+1.317(01) 


[0.0085,..., 0.0298] 


+0.0195(02) 


7m 


0.0085 


-0.0353(01) 


0.0285(01) 


-1.268(01) 


[0.0085,..., 0.0298] 


-0.0180(02) 


7p 


0.0085 


. n no s 1 / r\ i \ 

+0.0361(01) 


0.0285(01) 


.1 o /" o / n 1 \ 

+ 1.268(01) 


Tn nnoc n nonol 

[0.0085, . . . ,0.0298] 


. n n 1 o 1 / n i \ 

+0.0181(01) 


8m 


0.0020 


-0.0363(01) 


0.0280(01) 


-1.499(01) 


[0.0085,..., 0.0298] 


-0.0194(01) 


8p 


0.0020 


+0.0363(01) 


0.0274(01) 


+ 1.498(01) 


[0.0085,..., 0.0298] 


+0.0183(02) 


3m 


0.0180 


-0.0160(02) 


0.0218(01) 


-0.601(06) 


[0.0060,..., 0.0298] 


-0.0160(02) 


3p 


0.0180 


+0.0163(02) 


0.0219(01) 


+0.610(06) 


[0.0060,..., 0.0298] 


+0.0162(02) 


2m 


0.0085 


-0.0209(02) 


0.0182(01) 


-1.085(03) 


[0.0085,..., 0.0298] 


-0.0213(02) 


2p 


0.0085 


+0.0191(02) 


0.0170(02) 


+1.046(06) 


[0.0085,..., 0.0298] 


+0.0191(02) 


4m 


0.0085 


-0.0146(02) 


0.0141(01) 


-0.923(04) 


[0.0060,..., 0.0298] 


-0.0146(02) 


4p 


0.0085 


+0.0151(02) 


0.0144(01) 


+0.940(07) 


[0.0060,..., 0.0298] 


0.0151(02) 








j3 = 


2.10 






5m 


0.0078 


-0.00821(11) 


0.0102(01) 


-0.700(07) 


[0.0048,..., 0.0293] 


-0.0082(01) 


5p 


0.0078 


0.00823(08) 


0.0102(01) 


+0.701(05) 


[0.0048,..., 0.0293] 


+0.0082(01) 


4m 


0.0064 


-0.00682(13) 


0.0084(01) 


-0.706(09) 


in progress 




4p 


0.0064 


+0.00685(12) 


0.0084(01) 


+0.708(09) 


in progress 




3m 


0.0046 


-0.00585(08) 


0.0066(01) 


-0.794(07) 


[0.0025,..., 0.0297] 


-0.0059(01) 


3p 


0.0046 


+0.00559(14) 


0.0064(01) 


+0.771(13) 


[0.0025,..., 0.0297] 


+0.0056(01) 


2m 


0.0030 


-0.00403(14) 


0.0044(01) 


-0.821(17) 


in progress 




2P 


0.0030 


+0.00421(13) 


0.0045(01) 


+0.843(15) 


in progress 





Table 1: Overview of produced ensembles at /3 = 1.90, 1.95 and 2.10 



The basic ingredient of the calculation, due to Eq.(1.5), is the lattice RC estimator, 



yNp/m 

z r = 



Zp P//m (M Sea ' N tV m gSea,Np/m. |^val,Np/m gVal,Np/mj.^_2. ^ 



where j labels the valence quark polar mass and the momenta p are such that a 2 p 2 ranges from 0.5 
to 2.5 and that a 2 p^ / p 2 < 0.28, with pW = Y,pPp- Our analysis goes through the following steps. 



4 



Quark Bilinear RCs with Nf — 4 dynamical flavours 



David Palao 





0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 
(aA/, v!l1 + aJU™ 1 ) (aA/, v!l1 + aJU™ 1 ) 



(a) 



(b) 



Figure 1: j3 = 1.95, ensemble lp/m: (a) Goldstone pole subtraction fit applied separately on p and m for 
estimators of Zp 1 at p 2 ~ 9.5GeV 2 ; (b) extrapolation of the 0-averaged RC estimators to the chiral limit of 
Z q , Za, Zj and Zy at p 2 ~ 1 1 .5 GeV 2 . The situation is similar for other ensembles. 



(1) Subtract from the RC estimator the 0(a 2 g 2 ) cutoff effects at the chiral point, known from 
Ref. flu]] (2) Build the 0{a) improved estimator Z^ = Z^M^ -{Mj hl1 ,df^\,p 2 ^) = 



(3) Extrapolate to the chiral limit value, first in the valence and then in the sea sector. 

(4) Evaluate the RCs at a given renormalization scale after taking care of the residual lattice arte- 
facts according to the methods "Ml" and "M2" (see Ref. []7|] and discussion below). 

In step (2) chiral fit Ansatze are inspired to the mass parameter dependence expected from 
continuum QCD and the Symanzik analysis of lattice artefacts. For the valence chiral extrapolation 
we considered as fit functions linear combinations of constant, MJ al , (MJ a1 ) 2 , M v j al cos(6J iil ), and 
(MJ al cos{Bf 1 )) 2 . For the sea chiral limit linear combinations of M* ea , (M* ea ) 2 and cos(d] al ) (M° ea ) 2 
were considered. In the case of Zp, just before step (2) we remove the Goldstone pole contribution, 
which, depending directly on the lattice pseudoscalar meson mass, happens to be somewhat dif- 
ferent for estimators corresponding to opposite Q val J m -values. As for step (4), in the first method 
("Ml"), after bringing the RC-estimators to a common renormalization scale ipy^ = 1/a 2 ), we 
remove the remaining 0{a 2 p 2 ) discretization errors by a linear fit in p 2 . Here the fit range is 
1.5 < a 2 p 2 < 2.2. The second method ("M2") consists in simply taking the value of the RCs es- 
timators at some high momentum point kept fixed in physical units at all /3's. Here we choose 
p 2 = 12.0 ±0.5 GeV 2 . The two approaches yield RC values differing only by cutoff effects. 

In Fig. [TJ, for the example of the ensemble lp/m of j8 = 1.95, we show the Goldstone pole 
removal and the residual valence quark mass dependence in the analysis of Zp (panel (a)) and the 
extrapolation to the valence chiral limit for Z q , Za, Zt and Zy (panel (b)). The fit Ansatz is a linear 
function of the valence quark polar mass. We checked that results do not change significantly by 
using more complicated fit functions (involving higher mass powers or 6 ). Fig. ^| shows, for 
j8 = 1.95 the extrapolation to the sea chiral limit of Zp, Z$ (panel (a)) and Za, Zy (panel (b)). The 
fit Ansatz is a linear function of (M sea ) 2 . More elaborated fit functions give compatible results. 

For the typical (and important) cases of Z A and Zp we show for j8 = 1.95 (in Fig. ||) and 
j6 = 2. 1 (in Fig. 0) the residual dependence on a 2 p 2 of RC-estimators (in the case of Zp brought 
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0.70 - 
0.65 - 
0.60 - 
0.55 — 



Z s 
Z P 



Nil 



0.75 



0.70 



Za 
Zv 



0.50 0.65 

0e+00 2e-04 4e-04 6e-04 8e-04 le-03 le-03 0e+00 2e-04 4e-04 6e-04 8e-04 le-03 le-03 



(Mr? 

(a) 



(b) 



Figure 2: j3 = 1 .95: extrapolation to the sea chiral limit for (a) Z P and Z s ; (b) Z A and Z v (at p 2 ~ 1 1 .5 GeV 2 ). 




0.60 
0.55 
0.50 
0.45 
0.40 
0.35 



0.00 0.50 1.00 1.50 2.00 2.50 
(apf 

(a) 



uncorr 

2 




0.00 0.50 1.00 1.50 2.00 2.50 
(b) 



Figure 3: /3 = 1.95; dependence of Za (left panel) and Zp (right panel) on (a 2 p 2 ). Uncorrected and one-loop 
corrected (with two choices of the gauge coupling, bare and plaquette-boosted) RC estimators are shown. 



to a common renormalization scale (l/a(/$)) via three-loop evolution). The nice quality of the 
linear fit leading to the "Ml" RC-values is visible, while results of the "M2" type are obtained 
from data at a p in the range (1.8-2.0) and (1.10-1.23) for j3 = 1.95 and j3 =2.1, respectively. In 
each plot three different RC-estimators are considered, which differ from each other in the way the 
(beneficial) subtraction of 0(a 2 g 2 ) lattice artefacts of step (1) is carried out. 
In Table | we gather our preliminary results for the RCs at two values of the gauge coupling, 
j8 = 1 .95 and j8 = 2. 10. We present results obtained from the two methods described above, namely 
"Ml" and "M2". Perturbative contibutions 0(a 2 g 2 ) have been subtracted using the coupling con- 
stant estimate = 6/j8. Results for the RCs whose anomalous dimension is non-zero are given at 
the scale 1 /a in the RI/MOM scheme. Following ref. [[]]], we take a - 1 1^3=1.95(2.1) = 2.5(3.2) GeV. 
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Figure 4: /3 =2.10; dependence of (left panel) and Zp (right panel) on (a 2 p 2 ). Uncorrected and one-loop 
corrected (with two choices of the gauge coupling, bare and plaquette-boosted) RC estimators are shown. 



Method 


Z A 


Zv 


Z P {\/a) 


Zs(V«) 


ZpjZ s 


Z T (l/a) 


Z q (l/a) 








P = 


= 1.95 








Ml 


0.746(05) 


0.614(03) 


0.426(06) 


0.609(08) 


0.700(08) 


0.734(04) 


0.752(05) 


M2 


0.738(01) 


0.639(02) 


0.483(02) 


0.684(01) 


0.706(03) 


0.734(01) 


0.769(01) 








P = 


= 2.10 








Ml 


0.783(07) 


0.683(13) 


0.493(10) 


0.669(08) 


0.737(14) 


0.775(11) 


0.786(13) 


M2 


0.777(05) 


0.680(05) 


0.515(06) 


0.696(08) 


0.740(08) 


0.771(05) 


0.794(07) 



Table 2: Our preliminary results for quark bilinear RCs at J3 = 1 .95 and j5 =2.10. 
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